rm(list = ls())
#setwd("")

library(mediation)
library(haven)

data2 <- read_dta("hierarchy_replication.dta")


fit <- lm(Fwi_v2stfisccap2 ~  wi_usa_median +  wi_v2xcl_prpty + 
          wi_lnpop_wdi+ wi_lngdppc+ wi_polity2+ wi_polity2_2 + wi_ny_gdp_totl_rt_zs+ wi_cwyrs+ wi_c2+ wi_c3 +coldwar+wi_compete, data=data2)
esample.n <- nobs(fit)
esample<-rownames(as.matrix(resid(fit)))
data3<-data2[esample,]

##Model 2
m <- lm(wi_v2xcl_prpty ~ wi_usa_median + 
          wi_lnpop_wdi+ wi_lngdppc+ wi_polity2+ wi_polity2_2 + wi_ny_gdp_totl_rt_zs+ wi_cwyrs+ wi_c2+ wi_c3 +coldwar+wi_compete, data=data3)
##Model 3 
y <- lm(Fwi_v2stfisccap2 ~  wi_usa_median +  wi_v2xcl_prpty + 
           wi_lnpop_wdi+ wi_lngdppc+ wi_polity2+ wi_polity2_2 + wi_ny_gdp_totl_rt_zs+ wi_cwyrs+ wi_c2+ wi_c3 +coldwar+wi_compete, data=data3)

  esample.n <- nobs(fit)

  #Without bootstrapping (Not shown/reported)
  
 m.out<-mediate(m, y, treat = "wi_usa_median",
                 mediator = "wi_v2xcl_prpty")
 summary(m.out)
 plot(m.out, xlim=c(-0.02,0.12))
 
 #With bootstrapping (Figure 2c)
 m.out2<-mediate(m, y, boot=TRUE, sims=1000, treat = "wi_usa_median",
                mediator = "wi_v2xcl_prpty")
 summary(m.out2)
 plot(m.out2, xlim=c(-0.02,0.12) , xlab="Estimated Effect")
 
 sens.cont <- medsens(m.out2, rho.by = 0.05)
 summary(sens.cont)
 ##Figure B.1
 #Left Panel
 plot(sens.cont, sens.par = "rho")
 #Right Panel
 plot(sens.cont, sens.par = "R2")
 

 mediate.diagram(m.out2)